/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

// $Id: AliAODMuonReplicator0.cxx 56492 2012-05-15 18:42:47Z pcrochet $

//
// Implementation of a branch replicator 
// to produce slim muon and dimuon aods.
//
// This replicator is in charge of replicating the tracks,vertices,dimuons
// (and vzero and tzero) branches of the standard AOD into muon AODs 
// (AliAOD.Muons.root and AliAOD.Dimuons.root)
// 
// The tracks are filtered so that only muon tracks (and only muon tracks
// that pass the trackCut if present) make it to the output aods
//
// The vertices are filtered so that only the primary vertices make it
// to the output aods.
//
// The dimuons are recreated here, according to the set of tracks
// that pass the trackCut (that set may be the same as the input set,
// but to be 100% safe, we always recreate the dimuons).
// 
// Author: L. Aphecetche (Subatech)


class AliESDv0;
class AliESDVertex;
class AliAODVertex;
class AliAODRecoDecay;

#include "AliAODDimuon.h"
#include "AliAODEvent.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODTZERO.h"
#include "AliAODTrack.h"
#include "AliAODVZERO.h"
#include "AliAnalysisCuts.h"
#include "TF1.h"
#include "AliExternalTrackParam.h"
#include "AliESDv0.h"
#include "AliAODv0.h"
#include "AliPIDResponse.h"
//#include "AliAOD3LH.h"
#include <iostream>
#include <cassert>
#include "AliESDtrack.h"
#include "TObjArray.h"

#include "AliAODRecoDecay.h"
#include "AliAODRecoDecayHF.h"
#include "AliAODRecoDecayHF2Prong.h"

#include <TFile.h>
#include <TDatabasePDG.h>
#include <TString.h>
#include <TList.h>
#include "AliLog.h"
#include "AliVEvent.h"
#include "AliVVertex.h"
#include "AliVTrack.h"
#include "AliVertexerTracks.h"
#include "AliKFVertex.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrackCuts.h"
#include "AliAODEvent.h"
#include "AliAnalysisFilter.h"
#include "AliAnalysisVertexingHF.h"

#include "AliAODMuonReplicator0.h"
#include "TH1.h"
#include "TCanvas.h"

using std::cout;
using std::endl;

ClassImp(AliAODMuonReplicator0)

//_____________________________________________________________________________
AliAODMuonReplicator0::AliAODMuonReplicator0(const char* name, const char* title,
					     AliAnalysisCuts* trackCut,
					     AliAnalysisCuts* vertexCut,
					     Int_t mcMode, 
					     Int_t nsigmaTrk1, Int_t partType1,
					     Int_t nsigmaTrk2, Int_t partType2)
:AliAODBranchReplicator(name,title), 
  fTrackCut(trackCut), fTracks(0x0), 
  fVertexCut(vertexCut), fVertices(0x0), 
  fDimuons(0x0),
  fVZERO(0x0),
  fTZERO(0x0),
  fList(0x0),
  fMCParticles(0x0),
  fMCHeader(0x0),
  fnSigmaTrk1(nsigmaTrk1),
  fnSigmaTrk2(nsigmaTrk2),
  fpartType1(partType1),
  fpartType2(partType2),
  fMCMode(mcMode),
  fAODMapSize(0),
  fAODMap(0),
  fVertexerTracks(0x0),
  fBzkG(0.),
  fV1(0x0),
  fCosMin(),
  fDCAtracksMin(),
  fRmax(),
  fRmin(),
  fDNmin(),
  fDPmin(),
  fLabelMap(),
  fParticleSelected(),
  fSecVtxWithKF(kFALSE)
{
  // default ctor
}

//_____________________________________________________________________________
AliAODMuonReplicator0::~AliAODMuonReplicator0()
{
  // dtor
  delete fTrackCut;
  delete fVertexCut;
  delete fList;
}

//_____________________________________________________________________________
void AliAODMuonReplicator0::SelectParticle(Int_t i)
{
  // taking the absolute values here, need to take care 
  // of negative daughter and mother
  // IDs when setting!
  
  if (!IsParticleSelected(TMath::Abs(i)))
  {
    fParticleSelected.Add(TMath::Abs(i),1);    
  }
}

//_____________________________________________________________________________
Bool_t AliAODMuonReplicator0::IsParticleSelected(Int_t i)  
{
  // taking the absolute values here, need to take 
  // care with negative daughter and mother
  // IDs when setting!
  return (fParticleSelected.GetValue(TMath::Abs(i))==1);
}


//_____________________________________________________________________________
void AliAODMuonReplicator0::CreateLabelMap(const AliAODEvent& source)
{  
  //
  // this should be called once all selections are done 
  //
  
  fLabelMap.Delete();
  
  TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
  
  Int_t i(0);
  Int_t j(0);
  
  TIter next(mcParticles);
  
  while ( next() ) 
  {
    if (IsParticleSelected(i))
    {
      fLabelMap.Add(i,j++);
    }
    ++i;
  }
}

//_____________________________________________________________________________
Int_t AliAODMuonReplicator0::GetNewLabel(Int_t i) 
{
  // Gets the label from the new created Map
  // Call CreatLabelMap before
  // otherwise only 0 returned
  return fLabelMap.GetValue(TMath::Abs(i));
}

//_____________________________________________________________________________
void AliAODMuonReplicator0::FilterMC(const AliAODEvent& source)
{
  // Filter MC information

  fMCHeader->Reset();
  fMCParticles->Clear("C");

  AliAODMCHeader* mcHeader(0x0);
  TClonesArray* mcParticles(0x0);
  
  fParticleSelected.Delete();
  
  if ( fMCMode>=2 && !fTracks->GetEntries() ) return;
  // for fMCMode==1 we only copy MC information for events where there's at least one muon track
    
  mcHeader = static_cast<AliAODMCHeader*>(source.FindListObject(AliAODMCHeader::StdBranchName()));
  
  if ( mcHeader ) 
    {
      *fMCHeader = *mcHeader;
    }
  
  mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
  
  if ( mcParticles && fMCMode>=2 )
    {
      // loop on (kept) muon tracks to find their ancestors
      TIter nextMT(fTracks);
      AliAODTrack* mt;
    
      while ( ( mt = static_cast<AliAODTrack*>(nextMT()) ) )
	{
	  Int_t label = mt->GetLabel();
      
	  while ( label >= 0 ) 
	    {
	      SelectParticle(label);
	      AliAODMCParticle* mother = static_cast<AliAODMCParticle*>(mcParticles->UncheckedAt(label));
	      if (!mother)
		{
		  AliError("Got a null mother ! Check that !");
		  label = -1;
		}
	      else
		{
		  label = mother->GetMother();
		}
	    }
	}
    
      CreateLabelMap(source);
    
      // Actual filtering and label remapping (shamelessly taken for the implementation of AliAODHandler::StoreMCParticles)
      TIter nextMC(mcParticles);
      AliAODMCParticle* p;
      Int_t nmc(0);
      Int_t nmcout(0);
    
      while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
	{
	  AliAODMCParticle c(*p);
      
	  if ( IsParticleSelected(nmc) )
	    {
	      // 
	      Int_t d0 =  p->GetDaughter(0);
	      Int_t d1 =  p->GetDaughter(1);
	      Int_t m =   p->GetMother();
        
	      // other than for the track labels, negative values mean
	      // no daughter/mother so preserve it
        
	      if(d0<0 && d1<0)
		{
		  // no first daughter -> no second daughter
		  // nothing to be done
		  // second condition not needed just for sanity check at the end
		  c.SetDaughter(0,d0);
		  c.SetDaughter(1,d1);
		} else if(d1 < 0 && d0 >= 0) 
		{
		  // Only one daughter
		  // second condition not needed just for sanity check at the end
		  if(IsParticleSelected(d0))
		    {
		      c.SetDaughter(0,GetNewLabel(d0));
		    } else 
		    {
		      c.SetDaughter(0,-1);
		    }
		  c.SetDaughter(1,d1);
		}
	      else if (d0 > 0 && d1 > 0 )
		{
		  // we have two or more daughters loop on the stack to see if they are
		  // selected
		  Int_t d0tmp = -1;
		  Int_t d1tmp = -1;
		  for (int id = d0; id<=d1;++id)
		    {
		      if (IsParticleSelected(id))
			{
			  if(d0tmp==-1)
			    {
			      // first time
			      d0tmp = GetNewLabel(id);
			      d1tmp = d0tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same 
			    }
			  else d1tmp = GetNewLabel(id);
			}
		    }
		  c.SetDaughter(0,d0tmp);
		  c.SetDaughter(1,d1tmp);
		} else 
		{
		  AliError(Form("Unxpected indices %d %d",d0,d1));
		}
        
	      if ( m < 0 )
		{
		  c.SetMother(m);
		} else 
		{
		  if (IsParticleSelected(m)) 
		    {
		      c.SetMother(GetNewLabel(m));              
		    }
		  else 
		    {
		      AliError(Form("PROBLEM Mother not selected %d", m));              
		    }
		}
        
	      new ((*fMCParticles)[nmcout++]) AliAODMCParticle(c);
	    }
      
	  ++nmc;        
	}      
    
      // now remap the tracks...
    
      TIter nextTrack(fTracks);
      AliAODTrack* t;
    
      while ( ( t = static_cast<AliAODTrack*>(nextTrack()) ) )
	{
	  t->SetLabel(GetNewLabel(t->GetLabel()));
	}
    
    }
  else if ( mcParticles ) 
    {
      // simple copy of input MC particles to ouput MC particles
      TIter nextMC(mcParticles);
      AliAODMCParticle* p;
      Int_t nmcout(0);
    
      while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
	{
	  new ((*fMCParticles)[nmcout++]) AliAODMCParticle(*p);
	}
    }
  
  AliDebug(1,Form("input mc %d output mc %d",
                  mcParticles ? mcParticles->GetEntries() : 0,
                  fMCParticles ? fMCParticles->GetEntries() : 0));
  
}

//_____________________________________________________________________________
TList* AliAODMuonReplicator0::GetList() const
{
  // return (and build if not already done) our internal list of managed objects
  if (!fList)
    {
      // fTracks = new TClonesArray("AliAODTrack",30);
      // fTracks->SetName("tracks");    

      //      fTracks = new TClonesArray("AliAODv0",30);
      fTracks = new TClonesArray("AliAODRecoDecayHF2Prong",30);
      // fTracks->SetName("SecondaryVertices");    
      
      //  fTracks = new TClonesArray("AliESDv0",30);
      fTracks->SetName("SecondaryVertices");    

    
      fVertices = new TClonesArray("AliAODVertex",2);
      fVertices->SetName("vertices");    
    
      //    fDimuons = new TClonesArray("AliAODDimuon",30);
      fDimuons = new TClonesArray("AliAODTrack",30);
      fDimuons->SetName("Nuclei");
    
      fVZERO = new TClonesArray("AliAODTrack",30);
      fVZERO->SetName("DaughterTracks");

      fTZERO = new AliAODTZERO;
    
      fList = new TList;
      fList->SetOwner(kTRUE);
    
      fList->Add(fTracks);
      fList->Add(fVertices);
      fList->Add(fDimuons);
      fList->Add(fVZERO);
      
      //fList->Add(fTZERO);
    
      if ( fMCMode > 0 )
	{
	  fMCHeader = new AliAODMCHeader;    
	  fMCParticles = new TClonesArray("AliAODMCParticle",1000);
	  fMCParticles->SetName(AliAODMCParticle::StdBranchName());
	  fList->Add(fMCHeader);
	  fList->Add(fMCParticles);
	}
    }
  return fList;
}

//_____________________________________________________________________________
void AliAODMuonReplicator0::ReplicateAndFilter(const AliAODEvent& source)
{
  // Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
  
  printf("Sto eseguendo il muon Replicator...\n");
  
  assert(fTracks!=0x0);
  
  *fTZERO = *(source.GetTZEROData());
  
  fTracks->Clear("C");			
    
  fDimuons->Clear("C");

  fVZERO->Clear("C");

  // AliVEvent *event = new AliVEvent((AliAODEvent&)source);

  TIter next(source.GetTracks());

  Int_t ntracks(0);
  Int_t ndimuons(0);

  Int_t input(0);
  Double_t xdummy,ydummy;

  AliAODRecoDecayHF2Prong *io2Prong  = 0;

  TObjArray *twoTrackArray    = new TObjArray(2);
  Double_t dispersion;

  AliAODVertex *vtx = source.GetPrimaryVertex();
  
  // A Set of very loose cut for a weak strange decay
  
  fCosMin       = 0.97;
  fDCAtracksMin = 1;
  fRmax         = 200.;
  fRmin         = 0.1;
  fDNmin        = 0.05;
  fDPmin        = 0.05;

  //----------------------------------------------------------
  Int_t nindices=0;
  UShort_t *indices = 0;
  const Int_t entries = source.GetNumberOfTracks();
  //  AliCentrality* cent;

  Double_t pos[3],cov[6];
  vtx->GetXYZ(pos);
  vtx->GetCovarianceMatrix(cov);
  fV1 = new AliESDVertex(pos,cov,100.,100,vtx->GetName());
  cout<<"fV1 pricipal loop: "<<fV1<<endl;
  
  if(entries<=0) return;
  indices = new UShort_t[entries];
  memset(indices,0,sizeof(UShort_t)*entries);
  fAODMapSize = 100000;
  fAODMap = new Int_t[fAODMapSize];
  memset(fAODMap,0,sizeof(Int_t)*fAODMapSize);
  //  cent=((AliAODEvent&)source)->GetCentrality();
  
  //-------------------------------------------------------------

  AliAODRecoDecayHF   *rd = 0;

  if(vtx->GetNContributors()<1) {
    
    // SPD vertex cut
    vtx =source.GetPrimaryVertexSPD();
    
    if(vtx->GetNContributors()<1) {
      Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
      return; // NO GOOD VERTEX, SKIP EVENT 
    }
  }
  
  Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
  xPrimaryVertex=vtx->GetX();
  yPrimaryVertex=vtx->GetY();
  zPrimaryVertex=vtx->GetZ();  
  
  Double_t PrimaryVertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};

  Double_t lMagneticField=source.GetMagneticField();
  fBzkG=source.GetMagneticField();
  fVertexerTracks=new AliVertexerTracks(fBzkG);


  
  TF1 *foPion=new TF1("foPion", "[0]*([1]*TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3]) - 1 - TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3])*TMath::Log([2] + 1/TMath::Power((x/0.13957), [4])))",0.01,20);
  foPion->SetParameters(4.1,8.98482806165147636e+00,1.54000000000000005e-05,2.30445734159456084e+00,2.25624744086878559e+00);
  
  Double_t TrackNumber = source.GetNumberOfTracks();
  Int_t    countPionBuoni=0;
  
  //Tracks arrays
  
  TArrayI Track0(TrackNumber);     //Pions                                                                          
  Int_t nTrack0=0;
  
  TArrayI Track1(TrackNumber);        //Helium3
  Int_t nTrack1=0;
  
  
  ULong_t  status;
  Bool_t   isTPC, isTOF;
  
  for(Int_t j=0; j<TrackNumber; j++){
    
    AliAODTrack *aodtrack=source.GetTrack(j);
    
    //----------------------------------------
    // To be implemented into a method IsTrackAccepted 
    
    status  = (ULong_t)aodtrack->GetStatus();
    if((status&AliESDtrack::kTPCrefit)==0) continue;// Require tpcrefit
    
    //    if(aodtrack->GetTPCNcls() < 60) continue; // TPC Ncluster cut to be fixed
    if(aodtrack->GetTPCNcls() < 80) continue; // TPC Ncluster cut to be fixed
    if(aodtrack->Chi2perNDF() > 5) continue;
    
    //-----------------------------------------
    
    Double_t mom = aodtrack->GetDetPid()->GetTPCmomentum();
        
    if(mom<0.150)continue;
    
    // Double_t nSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType1)); //to be changed
    
    Double_t nSigmaNegPion = TMath::Abs((aodtrack->GetTPCsignal() - foPion->Eval(mom))/foPion->Eval(mom))/0.07;
    
    if ( (nSigmaNegPion <= fnSigmaTrk1) && aodtrack->Charge()==-1){
      Track0[nTrack0++]=j;
    }
    
    // Double_t nSigmaNuclei = fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType2);  //To be changed
    // if(nSigmaNuclei>-fnSigmaTrk2 && aodtrack->GetTPCsignal()<1000 && mom>0.2){ 
    
    Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
    
    if(aodtrack->GetTPCsignal() > triggerDeDx && aodtrack->GetTPCsignal()<5000 && mom>0.2 && aodtrack->Charge()==1){
      Track1[nTrack1++]=j;
      new((*fDimuons)[ndimuons++]) AliAODTrack(*aodtrack);
    }
  } 
  
  Track0.Set(nTrack0);
  Track1.Set(nTrack1);
  
  
  AliAODTrack *postrackAOD = 0;
  AliAODTrack *negtrackAOD = 0;
 
  AliESDtrack *postrack = 0;
  AliESDtrack *negtrack = 0;

  for (Int_t i=0; i<nTrack1; i++){                            //! He Tracks Loop

    Int_t Track1idx=Track1[i];
    
    postrackAOD=source.GetTrack(Track1idx);
     
    for (Int_t k=0; k <nTrack0 ; k++) {                           //! Pion Tracks Loop
      
      Int_t Track0idx=Track0[k];
      negtrackAOD=source.GetTrack(Track0idx);
      
      negtrack = new AliESDtrack((AliVTrack*)negtrackAOD);
      postrack = new AliESDtrack((AliVTrack*)postrackAOD);
      
      twoTrackArray->AddAt(negtrack,0);
      twoTrackArray->AddAt(postrack,1);
      
      Double_t dcap1n1 = postrack->GetDCA(negtrack,fBzkG,xdummy,ydummy);

      Double_t dcap1toPV = TMath::Abs(postrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));
      Double_t dcan1toPV = TMath::Abs(negtrack->GetD(xPrimaryVertex, yPrimaryVertex,fBzkG));

      if(dcap1n1>fDCAtracksMin)continue;
      if((xdummy+ydummy)>fRmax )continue;
      if((xdummy+ydummy)< fRmin)continue;
      
      if ( dcan1toPV < fDNmin)               
        if ( dcap1toPV < fDPmin) continue;   
      
     cout<<"dcap1n1: "<<dcap1n1<<endl;
 
      Int_t charge[2];
      charge[0]=negtrack->Charge();
      charge[1]=postrack->Charge();
    
      cout<<"Charge neg: "<<negtrack->Charge()<<endl;
      cout<<"Charge pos: "<<postrack->Charge()<<endl;
      
      AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray,dispersion,kTRUE);
      
      if(!vertexp1n1) {
      	// twoTrackArray->Clear();
      	// negtrack=0; 
      	continue; 
      }

      cout<<"vertexp1n1:  "<<vertexp1n1<<endl;

      io2Prong = Make2Prong(twoTrackArray,source,vertexp1n1,dcap1n1);
  
      if(io2Prong->CosPointingAngle()<fCosMin)continue;

      cout<<"pointing angle "<<io2Prong->CosPointingAngle()<<endl;
      
      rd =  new((*fTracks)[ntracks++]) AliAODRecoDecayHF2Prong(*io2Prong);
      rd->SetSecondaryVtx(vertexp1n1);
      vertexp1n1->SetParent(rd);
      //      AddRefs(vertexp1n1,rd,source,twoTrackArray);
      //vertexp1n1->Delete();  delete vertexp1n1;
	
    }
  }

  // twoTrackArray->Delete();  
  // delete twoTrackArray;
      
  //------------------------------
  
  // to be understand : do we need the information of the daughters? if yes we sould do someting like
  // $ALICE_ROOT/ANALYSIS/AliAnalysisTaskESDfilter.cxx :: Convert v0(AliESDevent) or not??
  // forse si se vogliamo poi fare dei controlli aggiuntivi sulle tracce
  // NB pero le abbiamo gia convertite, quindi sono solo piu da aggiungere in qualche modo
  
  //---------------------------------------------------------------------------------
  
  assert(fVertices!=0x0);
  fVertices->Clear("C");
  TIter nextV(source.GetVertices());
  AliAODVertex* v;
  Int_t nvertices(0);
  
  while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
    {
      if ( !fVertexCut || fVertexCut->IsSelected(v) ) 
	{
	  AliAODVertex* tmp = v->CloneWithoutRefs();
	  AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
	  
	  // to insure the main vertex retains the ncontributors information
	  // (which is otherwise computed dynamically from
	  // references to tracks, which we do not keep in muon aods...)
	  // we set it here
	  
	  copiedVertex->SetNContributors(v->GetNContributors()); 
	  delete tmp;
	}
    }
  
  printf("....Ho eseguito il muon Replicator...\n");
  
  AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d ndimuons=%d",
                  input,fTracks->GetEntries(),fVertices->GetEntries(),fDimuons->GetEntries())); 
  
  cout<<"....Sotto debug..."<<endl;
  
  // Finally, deal with MC information, if needed
  
  if ( fMCMode > 0 )
    {
      FilterMC(source);      
    }
  

  cout<<"Delete..."<<endl;
  cout<<"Delete 1"<<  endl;
  if(io2Prong) {delete io2Prong; io2Prong=NULL;}
  cout<<"Delete 2"<<  endl;
  twoTrackArray->Delete();  delete twoTrackArray;
  cout<<"Delete 3"<<  endl;
  // vtx->Delete();  delete vtx;
  cout<<"Delete 4"<<  endl;
  if(fV1) { delete fV1; fV1=NULL; }
  cout<<"Delete 5"<<  endl;
  if(fAODMap) { delete [] fAODMap; fAODMap=NULL; }
  cout<<"Delete 6"<<  endl;
  delete fVertexerTracks;
 // cout<<"Delete 7"<<  endl;
 //  postrackAOD->Delete();  delete  postrackAOD;
 //   cout<<"Delete 8"<<  endl;
 //  negtrackAOD->Delete();  delete  negtrackAOD;
 //  cout<<"Delete 9"<<  endl;
 //  postrack->Delete();  delete  postrack;
 //  cout<<"Delete 10"<<  endl;
 //  negtrack->Delete();  delete  negtrack;
  
  cout<<"Fine Replicator"<<endl;
}



//-----------------------------------------------------------------------------

AliAODVertex *AliAODMuonReplicator0::ReconstructSecondaryVertex(TObjArray *trkArray,
								Double_t &dispersion,Bool_t useTRefArray) const
{
  // Secondary vertex reconstruction with AliVertexerTracks or AliKFParticle
  //AliCodeTimerAuto("",0);
  cout<<"1"<<endl;
  AliESDVertex *vertexESD = 0;
  AliAODVertex *vertexAOD = 0;
  cout<<"2"<<endl;
  if(!fSecVtxWithKF) { // AliVertexerTracks
    cout<<"3"<<endl;
    fVertexerTracks->SetVtxStart(fV1);
    vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
    cout<<"4"<<endl;
    if(!vertexESD) return vertexAOD;
    cout<<"5"<<endl;
    // if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
    //   //AliDebug(2,"vertexing failed"); 
    //   delete vertexESD; vertexESD=NULL;
    //   return vertexAOD;
    // }
    cout<<"6"<<endl;
    Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
    if(vertRadius2>200.){
      // vertex outside beam pipe, reject candidate to avoid propagation through material
      delete vertexESD; vertexESD=NULL;
      return vertexAOD;
    }
    cout<<"7"<<endl;
  } else { // Kalman Filter vertexer (AliKFParticle)

    AliKFParticle::SetField(fBzkG);

    AliKFVertex vertexKF;

    Int_t nTrks = trkArray->GetEntriesFast();
    for(Int_t i=0; i<nTrks; i++) {
      AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
      AliKFParticle daughterKF(*esdTrack,211);
      vertexKF.AddDaughter(daughterKF);
    }
    vertexESD = new AliESDVertex(vertexKF.Parameters(),
				 vertexKF.CovarianceMatrix(),
				 vertexKF.GetChi2(),
				 vertexKF.GetNContributors());

  }
  cout<<"8"<<endl;
  // convert to AliAODVertex
  Double_t pos[3],cov[6],chi2perNDF;
  vertexESD->GetXYZ(pos); // position
  vertexESD->GetCovMatrix(cov); //covariance matrix
  chi2perNDF = vertexESD->GetChi2toNDF();
  dispersion = vertexESD->GetDispersion();
  delete vertexESD; vertexESD=NULL;
  cout<<"9"<<endl;
  Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);

  cout<<"------------------> Reconstruct vertexAOD: "<<vertexAOD<<endl;

  return vertexAOD;












  // //------------------------------------------------------------------------------

  // AliESDVertex *vertexESD = 0;
  // AliAODVertex *vertexAOD = 0;
 
  // if(!fSecVtxWithKF) { // AliVertexerTracks
  //   cout<<"fV1 inside loop: "<<fV1<<endl;
  //   cout<<"fVertexerTracks: "<<fVertexerTracks<<endl;
  //   fVertexerTracks->SetVtxStart(fV1);
  //   vertexESD = (AliESDVertex*)fVertexerTracks->VertexForSelectedESDTracks(trkArray);
  //   cout<<vertexESD->GetXv()<<" "<<vertexESD->GetYv()<<" "<<vertexESD->GetZv()<<endl;
  //   // if(!vertexESD) 
  //   //   return vertexAOD;
  //   cout<<"6"<<endl;
    
  //   if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
  //     cout<<"7"<<endl;
  //     //AliDebug(2,"vertexing failed"); 
  //     //      delete vertexESD; vertexESD=NULL;
  //     cout<<"8"<<endl;
  //     cout<<"vertexAOD: "<<vertexAOD<<endl;
  //     return vertexAOD;
  //   }

  //   cout<<vertexESD->GetXv()<<" "<<vertexESD->GetYv()<<" "<<vertexESD->GetZv()<<endl;
  //   Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
  //   cout<<"Radius: "<<vertRadius2<<endl;
  //   if(vertRadius2>200.){
  //     // vertex outside beam pipe, reject candidate to avoid propagation through material
  //     delete vertexESD; vertexESD=NULL;
  //     return vertexAOD;
  //   }
  //   cout<<"9"<<endl;
  // } else { // Kalman Filter vertexer (AliKFParticle)

  //   AliKFParticle::SetField(fBzkG);

  //   AliKFVertex vertexKF;

  //   Int_t nTrks = trkArray->GetEntriesFast();
  //   for(Int_t i=0; i<nTrks; i++) {
  //     AliESDtrack *esdTrack = (AliESDtrack*)trkArray->At(i);
  //     AliKFParticle daughterKF(*esdTrack,211);
  //     vertexKF.AddDaughter(daughterKF);
  //   }
  //   vertexESD = new AliESDVertex(vertexKF.Parameters(),
  // 				 vertexKF.CovarianceMatrix(),
  // 				 vertexKF.GetChi2(),
  // 				 vertexKF.GetNContributors());

  // }
  
  // // convert to AliAODVertex
  // Double_t pos[3],cov[6],chi2perNDF;
  // vertexESD->GetXYZ(pos); // position
  // vertexESD->GetCovMatrix(cov); //covariance matrix
  // chi2perNDF = vertexESD->GetChi2toNDF();
  // dispersion = vertexESD->GetDispersion();
  // delete vertexESD; vertexESD=NULL;
  
  // Int_t nprongs= (useTRefArray ? 0 : trkArray->GetEntriesFast());
  // vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
  // return vertexAOD;

}

//-----------------------------------------------------------------------------
// AliAODRecoDecayHF2Prong *AliAODMuonReplicator0::Make2Prong(
// 							   TObjArray *twoTrackArray,AliVEvent &event,
// 							   AliAODVertex *secVert,Double_t dca,Int_t *charge[2]) 
// AliAODRecoDecayHF2Prong *AliAODMuonReplicator0::Make2Prong(
// 							   TObjArray *twoTrackArray,AliVEvent &event,
// 							   AliAODVertex *secVert,Double_t dca) 

// AliAODRecoDecayHF2Prong *AliAODMuonReplicator0::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
// 							   AliAODVertex *secVert,Double_t dcap1n1)

AliAODRecoDecayHF2Prong* AliAODMuonReplicator0::Make2Prong(TObjArray *twoTrackArray,const AliAODEvent &evento,
							   AliAODVertex *secVert,Double_t dca)

{ 

  cout<<"Inside make 2 prong"<<endl;

  Int_t charge[2];
  charge[0]=-1;
  charge[1]=2;
      
  // From AliAnalysisVertexingHF.cxx Method:Make2Prongs
  
  cout<< charge[0]<<" "<< charge[1]<<endl;
      
  fBzkG = evento.GetMagneticField();
      
  cout<<"fBzkG: "<<fBzkG<<endl;
 cout<<"twoTrackArray: "<<twoTrackArray<<endl;
  Double_t px[2],py[2],pz[2],d0[2],d0err[2];
  AliESDtrack *postrack = (AliESDtrack*)twoTrackArray->UncheckedAt(0);
  AliESDtrack *negtrack = (AliESDtrack*)twoTrackArray->UncheckedAt(1);
  
  cout<<"negtrack: "<<negtrack<<" postrack: "<<postrack<<endl;
  cout<<"kVeryBig: "<<kVeryBig<<endl;
  cout<<"secVert: "<<secVert<<endl;

  // // propagate tracks to secondary vertex, to compute inv. mass
  postrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
  negtrack->PropagateToDCA(secVert,fBzkG,kVeryBig);
  
  Double_t momentum[3];
  postrack->GetPxPyPz(momentum);
  px[0] = charge[1]*momentum[0]; py[0] = charge[1]*momentum[1]; pz[0] = charge[1]*momentum[2]; 
  cout<< px[0] <<" "<< " "<< py[0] << " "<< pz[0]<<endl;
  //  px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2]; 
  negtrack->GetPxPyPz(momentum);
  px[1] = charge[0]*momentum[0]; py[1] = charge[0]*momentum[1]; pz[1] = charge[0]*momentum[2]; 
  cout<< px[1] <<" "<< " "<< py[1] << " "<< pz[1]<<endl;
  //px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2]; 
  

  // primary vertex to be used by this candidate
  AliAODVertex *primVertexAOD  = evento.GetPrimaryVertex();
  cout<<"primVertexAOD "<<primVertexAOD<<endl;
  if(!primVertexAOD) return 0x0;
      
  Double_t d0z0[2],covd0z0[3];
  //  postrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
  d0[0] = d0z0[0];
  d0err[0] = TMath::Sqrt(covd0z0[0]);
  // negtrack->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
  d0[1] = d0z0[0];
  d0err[1] = TMath::Sqrt(covd0z0[0]);
  
  // create the object AliAODRecoDecayHF2Prong
  //  AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dcap1n1);
  AliAODRecoDecayHF2Prong *the2Prong = new AliAODRecoDecayHF2Prong(secVert,px,py,pz,d0,d0err,dca);
  cout<<"the2Prong "<<the2Prong<<endl;
  the2Prong->SetOwnPrimaryVtx(primVertexAOD);
  
  UShort_t id[2]={(UShort_t)postrack->GetID(),(UShort_t)negtrack->GetID()};
  the2Prong->SetProngIDs(2,id);
  //delete primVertexAOD; primVertexAOD=NULL;
 
  if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar 
      
    AddDaughterRefs(secVert,(AliAODEvent&)evento,twoTrackArray);
      
  }
  
  return the2Prong;  
}

//----------------------------------------------------------------------------
void AliAODMuonReplicator0::AddDaughterRefs(AliAODVertex *v,
					    const AliAODEvent &event,
					    const TObjArray *trkArray) const
{
  // Add the AOD tracks as daughters of the vertex (TRef)
  //AliCodeTimerAuto("",0);
  cout<<"Inside  AddDaughterRefs "<<endl;

  Int_t nDg = v->GetNDaughters();
  
  cout<<"nDg "<<nDg<<endl;

  TObject *dg = 0;
  if(nDg) dg = v->GetDaughter(0);
  cout<<"1 "<<1<<endl;
  if(dg) return; // daughters already added
  cout<<"2 "<<2<<endl;
  Int_t nTrks = trkArray->GetEntriesFast();
  cout<<"nTrks "<<nTrks<<endl;
  AliExternalTrackParam *track = 0;
  AliAODTrack *aodTrack = 0;
  Int_t id;
  
  for(Int_t i=0; i<nTrks; i++) {
    track = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
    cout<<"track inside loop: "<<track<<endl;
    id = (Int_t)track->GetID();
    printf("---> %d\n",id);
    if(id<0) continue; // this track is a AliAODRecoDecay
    cout<<aodTrack<<endl;
    cout<<"1 "<<fAODMap[id]<<endl;
    aodTrack = (AliAODTrack*)event.GetTrack(fAODMap[id]);
    cout<<"2 "<<fAODMap[id]<<endl;
    v->AddDaughter(aodTrack);
    cout<<"End Loop"<<endl;
  }
  cout<<"END AddDaughterRefs"<<endl;
  return;
}
//----------------------------------------------------------------------------
	
void AliAODMuonReplicator0::AddRefs(AliAODVertex *v,AliAODRecoDecayHF *rd,
				    const AliAODEvent &event,
				    const TObjArray *trkArray) const
{
  // Add the AOD tracks as daughters of the vertex (TRef)
  // Also add the references to the primary vertex and to the cuts
  //AliCodeTimerAuto("",0);
  
  AddDaughterRefs(v,event,trkArray);
  rd->SetPrimaryVtxRef((AliAODVertex*)event.GetPrimaryVertex());
  return;
}	

//-----------------------------------------------------------------------------

// AliAODVertex* AliAODMuonReplicator0::PrimaryVertex(const TObjArray *trkArray,
// 						   AliAODEvent &event) const
// {
//   // Returns primary vertex to be used for this candidate
//   //AliCodeTimerAuto("",0);

//   AliESDVertex *vertexESD = 0;
//   AliAODVertex *vertexAOD = 0;


//   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) { 
//     // primary vertex from the input event
    
//     vertexESD = new AliESDVertex(*fV1);

//   } else {
//     // primary vertex specific to this candidate

//     Int_t nTrks = trkArray->GetEntriesFast();
//     AliVertexerTracks *vertexer = new AliVertexerTracks(event->GetMagneticField());

//     if(fRecoPrimVtxSkippingTrks) { 
//       // recalculating the vertex
      
//       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
// 	Float_t diamondcovxy[3];
// 	event->GetDiamondCovXY(diamondcovxy);
// 	Double_t pos[3]={event->GetDiamondX(),event->GetDiamondY(),0.};
// 	Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
// 	AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
// 	vertexer->SetVtxStart(diamond);
// 	delete diamond; diamond=NULL;
// 	if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter")) 
// 	  vertexer->SetOnlyFitter();
//       }
//       Int_t skipped[1000];
//       Int_t nTrksToSkip=0,id;
//       AliExternalTrackParam *t = 0;
//       for(Int_t i=0; i<nTrks; i++) {
// 	t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
// 	id = (Int_t)t->GetID();
// 	if(id<0) continue;
// 	skipped[nTrksToSkip++] = id;
//       }
//       // TEMPORARY FIX
//       // For AOD, skip also tracks without covariance matrix
//       if(fInputAOD) {
// 	Double_t covtest[21];
// 	for(Int_t j=0; j<event->GetNumberOfTracks(); j++) {
// 	  AliVTrack *vtrack = (AliVTrack*)event->GetTrack(j);
// 	  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
// 	    id = (Int_t)vtrack->GetID();
// 	    if(id<0) continue;
// 	    skipped[nTrksToSkip++] = id;
// 	  }
// 	}
//       }
//       for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
//       //
//       vertexer->SetSkipTracks(nTrksToSkip,skipped);
//       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
      
//     } else if(fRmTrksFromPrimVtx && nTrks>0) { 
//       // removing the prongs tracks
      
//       TObjArray rmArray(nTrks);
//       UShort_t *rmId = new UShort_t[nTrks];
//       AliESDtrack *esdTrack = 0;
//       AliESDtrack *t = 0;
//       for(Int_t i=0; i<nTrks; i++) {
// 	t = (AliESDtrack*)trkArray->UncheckedAt(i);
// 	esdTrack = new AliESDtrack(*t);
// 	rmArray.AddLast(esdTrack);
// 	if(esdTrack->GetID()>=0) {
// 	  rmId[i]=(UShort_t)esdTrack->GetID();
// 	} else {
// 	  rmId[i]=9999;
// 	}
//       }
//       Float_t diamondxy[2]={event->GetDiamondX(),event->GetDiamondY()};
//       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
//       delete [] rmId; rmId=NULL;
//       rmArray.Delete();
      
//     }

//     if(!vertexESD) return vertexAOD;
//     if(vertexESD->GetNContributors()<=0) { 
//       //AliDebug(2,"vertexing failed"); 
//       delete vertexESD; vertexESD=NULL;
//       return vertexAOD;
//     }

//     delete vertexer; vertexer=NULL;

//   }

//   // convert to AliAODVertex
//   Double_t pos[3],cov[6],chi2perNDF;
//   vertexESD->GetXYZ(pos); // position
//   vertexESD->GetCovMatrix(cov); //covariance matrix
//   chi2perNDF = vertexESD->GetChi2toNDF();
//   delete vertexESD; vertexESD=NULL;

//   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);

//   return vertexAOD;
// }


// //---------------------------------------------------------------------------
//______________________________________________________________________________
TClonesArray& AliAODMuonReplicator0::fV0s()
{
  TClonesArray pippo("AliAODv0",30);//= new TClonesArray("AliAODv0",30);
  return pippo;
  //  return *(AODEvent()->GetV0s());
  //  return *(source->GetV0s());
  // return *0x0;
}

void AliAODMuonReplicator0::Terminate(){

}
